/*

Replication file for 
Kurt Schmidheiny and Sebastian Siegloch, 
"On Event Studies and Distributed-Lags in Two-Way Fixed Effects Models: 
Identification, Equivalence, and Generalization", 
Journal of Applied Econometrics, 2023.

(c) Kurt Schmidheiny and Sebastian Siegloch (2020, 2023)

Version 1 (May 27, 2020)
Version 2 (Jan 30, 2023)

*/
	
version 17

* clear workspace
clear
estimates clear

* uncomment to install package to plot estimated coefficients
// ssc install coefplot, replace

* set local directory
* cd "[yourpath]"

* open data from Baker and Fradkin (2017, henceforth BF2017)
use bf2017.dta

* create period variable from year and month in Stata format (Jan 1960 = 0)
generate yearmonth = 12*(year-1960)+month-1
format yearmonth %tm
label var yearmonth "Serial number for month and year (Jan 1960 = 0)"

* log of dependent variable
generate ln_GJSI = ln(GJSI)
label var ln_GJSI "log of GJSI"

* covariate used in BF2017
generate frac_total_ui = (cont_claims+init_claims)/population
label var frac_total_ui "Unemployment insurance claims in state and month divided by state population"

order state state_abb year month yearmonth

* set panel
xtset state yearmonth

* ------------------------------------------------------------------------------
* Figure B.1, Panel A, left
* multiple treatments of identical intensities
* replicates Baker and Fradkin, Table 8, col (3)
* ------------------------------------------------------------------------------	

* Baker and Fradkin use whole sample of dependent variable
* they assume that there was no treatment before and after observed sample

* create treatment adoption indicator
* increase by at least 13 weeks
generate D_PBD_incr = (D.PBD>=13) if !missing(D.PBD)
sum D_PBD_incr
* replace missing first treatment with zero (implicitly assumed by BF2017)
replace D_PBD_incr = 0 if missing(D_PBD_incr)

* create 3 leads and 4 lags
foreach lag in F3 F2 F L L2 L3 L4 {
	generate `lag'_D_PBD_incr = `lag'.D_PBD_incr
	* replace missing lags with zero (implicitly assumed by BF2017)
	replace `lag'_D_PBD_incr = 0 if missing(`lag'_D_PBD_incr)
}

* view resulting design matrix
// browse state year month yearmonth D_PBD_* *_D_PBD_*

* estimate with fixed effects
xtreg ln_GJSI ///
	F3_D_PBD_incr F2_D_PBD_incr F_D_PBD_incr ///
	D_PBD_incr ///
	L_D_PBD_incr L2_D_PBD_incr L3_D_PBD_incr L4_D_PBD_incr ///
	frac_total_ui i.yearmonth ///
	if year <= 2011 ///
	, fe vce(cluster state)

estimates store beta_incr_bf
local obs=e(N)
local T=e(g_avg)
local N=e(N_g)
sum year if e(sample)
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & e(sample)
local month_start = r(min)
sum month if year==`year_end' & e(sample)
local month_end = r(max)

* plot coefficients
coefplot (beta_incr_bf, ciopts(recast(rcap)) recast(connected) level(95)) ///
  	, ///
  	vertical ///
  	keep (*D_PBD_incr) ///
  	coeflabel(F3_D_PBD_incr = -3 F2_D_PBD_incr = -2 F_D_PBD_incr = -1 D_PBD_incr = 0 L_D_PBD_incr = 1 L2_D_PBD_incr = 2 L3_D_PBD_incr = 3 L4_D_PBD_incr = 4) ///
  	xline(3.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medlarge)) ///
 	yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medlarge)) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medlarge)) ytitle ("Effect on log search intensity" " ", size(medlarge)) ///
	legend(off) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medlarge))

graph export "Fig_B1_A_left.pdf", replace

* ------------------------------------------------------------------------------
* Figure B.1, Panel A, right
* multiple treatments of identical intensities
* replicates Baker and Fradkin, Table 8, col (5)
* ------------------------------------------------------------------------------	

* create treatment adoption indicator
* decrease by at least 7 weeks
generate D_PBD_decr = (D.PBD<=-7) if !missing(D.PBD)

* create 3 leads and 4 lags
foreach lag in F3 F2 F L L2 L3 L4 {
	generate `lag'_D_PBD_decr = `lag'.D_PBD_decr
}

* estimate
xtreg ln_GJSI ///
	F3_D_PBD_decr F2_D_PBD_decr F_D_PBD_decr ///
	D_PBD_decr ///
	L_D_PBD_decr L2_D_PBD_decr L3_D_PBD_decr L4_D_PBD_decr ///
	frac_total_ui i.yearmonth ///
	if year >= 2012 ///
	, fe vce(cluster state)

estimates store beta_decr_bf
local obs=e(N)
local T=e(g_avg)
local N=e(N_g)
sum year if e(sample)
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & e(sample)
local month_start = r(min)
sum month if year==`year_end' & e(sample)
local month_end = r(max)

* plot coefficients
coefplot (beta_decr_bf, ciopts(recast(rcap)) recast(connected) level(95)) ///
  	, ///
  	vertical ///
  	keep (*D_PBD_decr) ///
  	coeflabel(F3_D_PBD_decr = -3 F2_D_PBD_decr = -2 F_D_PBD_decr = -1 D_PBD_decr = 0 L_D_PBD_decr = 1 L2_D_PBD_decr = 2 L3_D_PBD_decr = 3 L4_D_PBD_decr = 4) ///
  	xline(3.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medlarge)) ///
 	yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medlarge)) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medlarge)) ytitle ("Effect on log search intensity" " ", size(medlarge)) ///
	legend(off) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medlarge))

graph export "Fig_B1_A_right.pdf", replace

* ------------------------------------------------------------------------------
* Figure B.1, Panel B, left
* multiple treatments of identical intensities (Schmidheiny/Siegloch case 2)
* effect window from -3 to +4
* estimation with distributed-lags in levels
* ------------------------------------------------------------------------------

* generate treatment status from treatment adoption indicator for increases
* cumulative sum calculated recursively, arbitrary starting value set to zero
generate PBD_incr = 0 if yearmonth==552 // first month, Jan 2006
replace  PBD_incr = L.PBD_incr+D_PBD_incr if yearmonth>552 
sum PBD_incr

* estimate distributed-lag model in levels with 2 leads and 4 lags
xtreg ln_GJSI ///
	F2.PBD_incr F.PBD_incr ///
	PBD_incr ///
	L.PBD_incr L2.PBD_incr L3.PBD_incr L4.PBD_incr ///
	frac_total_ui i.yearmonth ///
	if year <= 2011 ///
	, fe vce(cluster state)

estimates store gamma_incr_dl_fe
local obs=e(N)
local T=e(g_avg)
local N=e(N_g)
sum year if e(sample)
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & e(sample)
local month_start = r(min)
sum month if year==`year_end' & e(sample)
local month_end = r(max)

* cumulative sum of gamma coefficients starting at zero in period -1
nlcom ///
	(-(_b[F1.PBD_incr]+_b[F2.PBD_incr])) ///
	(-_b[F1.PBD_incr]) ///
	(0) ///
	(_b[PBD_incr]) ///
	(_b[PBD_incr]+_b[L1.PBD_incr]) ///
	(_b[PBD_incr]+_b[L1.PBD_incr]+_b[L2.PBD_incr]) ///
	(_b[PBD_incr]+_b[L1.PBD_incr]+_b[L2.PBD_incr]+_b[L3.PBD_incr]) ///
	(_b[PBD_incr]+_b[L1.PBD_incr]+_b[L2.PBD_incr]+_b[L3.PBD_incr]+_b[L4.PBD_incr]) ///
	, post level(95)

estimates store beta_incr_dl_fe

* plot beta coefficients
coefplot (beta_incr_dl_fe, ciopts(recast(rcap)) recast(connected) level(95)) ///
  	, ///
  	vertical ///
  	keep (_nl_*) ///
  	coeflabel(_nl_1 = -3 _nl_2 = -2 _nl_3 = -1 _nl_4 = 0 _nl_5 = 1 _nl_6 = 2 _nl_7 = 3 _nl_8 = 4) ///
  	xline(3.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medlarge)) ///
 	yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medlarge)) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medlarge)) ytitle ("Effect on log search intensity" " ", size(medlarge)) ///
	legend(off) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medlarge))
	
graph export "Fig_B1_B_left.pdf", replace

* ------------------------------------------------------------------------------
* Figure B.1, Panel B, right
* multiple treatments of identical intensities (Schmidheiny/Siegloch case 2)
* effect window from -3 to +4
* estimation with distributed-lags in levels
* ------------------------------------------------------------------------------

* generate treatment status from treatment adoption indicator for decreases
* cumulative sum calculated recursively, arbitrary starting value set to zero
generate PBD_decr = 0 if yearmonth==552 // first month, Jan 2006
replace  PBD_decr = L.PBD_decr+D_PBD_decr if yearmonth>552
sum PBD_decr

* estimate distributed-lag model in levels with 2 leads and 4 lags
xtreg ln_GJSI ///
	F2.PBD_decr F.PBD_decr ///
	PBD_decr ///
	L.PBD_decr L2.PBD_decr L3.PBD_decr L4.PBD_decr ///
	frac_total_ui i.yearmonth ///
	if year >= 2012 ///
	, fe vce(cluster state)

local obs=e(N)
local T=e(g_avg)
local N=e(N_g)
sum year if e(sample)
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & e(sample)
local month_start = r(min)
sum month if year==`year_end' & e(sample)
local month_end = r(max)

* cumulative sum of gamma coefficients starting at zero in period -1
nlcom ///
	(-(_b[F1.PBD_decr]+_b[F2.PBD_decr])) ///
	(-_b[F1.PBD_decr]) ///
	(0) ///
	(_b[PBD_decr]) ///
	(_b[PBD_decr]+_b[L1.PBD_decr]) ///
	(_b[PBD_decr]+_b[L1.PBD_decr]+_b[L2.PBD_decr]) ///
	(_b[PBD_decr]+_b[L1.PBD_decr]+_b[L2.PBD_decr]+_b[L3.PBD_decr]) ///
	(_b[PBD_decr]+_b[L1.PBD_decr]+_b[L2.PBD_decr]+_b[L3.PBD_decr]+_b[L4.PBD_decr]) ///
	, post level(95)

estimates store beta_decr_dl_fe

* plot beta coefficients
coefplot (beta_decr_dl_fe, ciopts(recast(rcap)) recast(connected) level(95)) ///
  	, ///
  	vertical ///
  	keep (_nl_*) ///
  	coeflabel(_nl_1 = -3 _nl_2 = -2 _nl_3 = -1 _nl_4 = 0 _nl_5 = 1 _nl_6 = 2 _nl_7 = 3 _nl_8 = 4) ///
  	xline(3.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medlarge)) ///
 	yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medlarge)) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medlarge)) ytitle ("Effect on log search intensity" " ", size(medlarge)) ///
	legend(off) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medlarge))
	
graph export "Fig_B1_B_right.pdf", replace

* ------------------------------------------------------------------------------
* Figure B.1, Panel B, left 
* alternative estimation with event study specification
* ------------------------------------------------------------------------------

* drop treatment adoption indicators from above with assumed zeros for NAs
drop *D_PBD_incr

* create treatment adoption indicator
* increase by at least 13 weeks
generate D_PBD_incr = (D.PBD>=13) if !missing(D.PBD)

* treatment adoption indicators with 3 leads and 4 lags
foreach lag in F3 F2 F L L2 L3 L4 {
	generate `lag'_D_PBD_incr = `lag'.D_PBD_incr
}

* The first observation of the treatment status in the data is 2006-01. The first
* treatment adoption indicator is observed in 2006-02. The first observation of
* the dependent variable that can be included in the estimation is 4-1=3 periods
* later, i.e. 2006-5 (Schmidheiny/Siegloch Remark 4). The unobserved 4th lag 
* of the treatment adopton indicator for 2006-5 can be set to an arbitrary value,
* e.g. zero. 
replace L4_D_PBD_incr = 0 if yearmonth==(556) // yearmonth 556 is 2006-05
* Note that in unbalanced panel, this operation may be unit specific.

* The last observation of the dependent variable in the crisis sample is 2011-12.
* There are more than 3 leads observed as the treatment adoption indicator is 
* observed up to 2015-12. No need to set arbitrary starting value.

* generate binned endpoint at lag 4 according to eq. (5)
* cumulative sum of lag 4 starting at 2006-05 (yearmonth 556)
sort state yearmonth
generate L4bin_D_PBD_incr = L4_D_PBD_incr if yearmonth==556
replace  L4bin_D_PBD_incr = L4bin_D_PBD_incr[_n-1] + L4_D_PBD_incr if yearmonth>556

* generate binned endpoint at lead 3 according to eq. (5)
* downward cumulative sum of lead 3 starting at 2011-12 (yearmonth 623)
gsort state -yearmonth
generate F3bin_D_PBD_incr = F3_D_PBD_incr if yearmonth==623
replace  F3bin_D_PBD_incr = F3bin_D_PBD_incr[_n-1] + F3_D_PBD_incr if yearmonth<623
sort state yearmonth

* view binned and not binned endpoints
* browse state yearmonth L4_D_PBD_incr L4bin_D_PBD_incr F3_D_PBD_incr F3bin_D_PBD_incr 

* estimate event-study model in levels with fixed effects
* with 4 lags, 3 leads, and normalization at -1
xtreg ln_GJSI ///
	F3bin_D_PBD_incr ///
	F2_D_PBD_incr ///
	D_PBD_incr ///
	L_D_PBD_incr L2_D_PBD_incr L3_D_PBD_incr ///
	L4bin_D_PBD_incr ///
	frac_total_ui i.yearmonth ///
	if year <= 2011 ///
	, fe vce(cluster state)

local obs=e(N)
local T=e(g_avg)
local N=e(N_g)
sum year if e(sample)
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & e(sample)
local month_start = r(min)
sum month if year==`year_end' & e(sample)
local month_end = r(max)

* add zero at period -1 to estimates
nlcom ///
	(_b[F3bin_D_PBD_incr]) ///
	(_b[F2_D_PBD_incr]) ///
	(0) ///
	(_b[D_PBD_incr]) ///
	(_b[L_D_PBD_incr]) ///
	(_b[L2_D_PBD_incr]) ///
	(_b[L3_D_PBD_incr]) ///
	(_b[L4bin_D_PBD_incr]) ///
	, post level(95)

estimates store beta_incr_es_fe

* plot beta coefficients
coefplot (beta_incr_es_fe, ciopts(recast(rcap)) recast(connected) level(95)) ///
  	, ///
  	vertical ///
  	keep (_nl_*) ///
  	coeflabel(_nl_1 = -3 _nl_2 = -2 _nl_3 = -1 _nl_4 = 0 _nl_5 = 1 _nl_6 = 2 _nl_7 = 3 _nl_8 = 4) ///
  	xline(3.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medlarge)) ///
 	yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medlarge)) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medlarge)) ytitle ("Effect on log search intensity" " ", size(medlarge)) ///
	legend(off) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medlarge))

* identical sample, point estimates and confidence bounds as in distributed-lag specification

* ------------------------------------------------------------------------------
* Figure B.2, left
* multiple treatments of varying intensities (Schmidheiny/Siegloch case 4)
* effect window from -3 to +4
* estimation with distributed-lags in levels
* full sample
* ------------------------------------------------------------------------------

* descsribe treatment status
sum PBD

* estimate distributed-lag model in levels with fixed effects with 2 leads and 4 lags
xtreg ln_GJSI ///
	F2.PBD F.PBD ///
	PBD ///
	L1.PBD L2.PBD L3.PBD L4.PBD ///
	frac_total_ui i.yearmonth ///
	, fe vce(cluster state)

estimates store gamma_dl_fe_full
local obs=e(N)
local T=e(g_avg)
local N=e(N_g)
sum year if e(sample)
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & e(sample)
local month_start = r(min)
sum month if year==`year_end' & e(sample)
local month_end = r(max)

* cumulative sum of gamma coefficients starting at zero in period -1
nlcom ///
	(-(_b[F1.PBD]+_b[F2.PBD])) ///
	(-_b[F1.PBD]) ///
	(0) ///
	(_b[PBD]) ///
	(_b[PBD]+_b[L1.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]+_b[L3.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]+_b[L3.PBD]+_b[L4.PBD]) ///
	, post level(95)

estimates store beta_dl_fe_full

* plot beta coefficients
coefplot (beta_dl_fe_full, ciopts(recast(rcap)) recast(connected) level(95)) ///
  	, ///
  	vertical ///
  	keep (_nl_*) ///
  	coeflabel(_nl_1 = -3 _nl_2 = -2 _nl_3 = -1 _nl_4 = 0 _nl_5 = 1 _nl_6 = 2 _nl_7 = 3 _nl_8 = 4) ///
  	xline(3.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medlarge)) ///
 	ylabel(-0.006(0.002)0.002) yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medlarge)) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medlarge)) ytitle ("Effect on log search intensity" " ", size(medlarge)) ///
	legend(off) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medlarge))

graph export "Fig_B2_left.pdf", replace
	
* ------------------------------------------------------------------------------
* Figure B.2, right
* also reported in Figure B.4, left (FE)
* multiple treatments of varying intensities (Schmidheiny/Siegloch case 4)
* effect window from -3 to +4
* estimation with distributed-lags in levels
* crisis sample
* ------------------------------------------------------------------------------

* descsribe treatment status
sum PBD

* estimate distributed-lag model in levels with fixed effects with 2 leads and 4 lags
xtreg ln_GJSI ///
	F2.PBD F.PBD ///
	PBD ///
	L1.PBD L2.PBD L3.PBD L4.PBD ///
	frac_total_ui i.yearmonth ///
	if year <= 2011 ///
	, fe vce(cluster state)

estimates store gamma_dl_fe
generate sample_dl_fe_crisis = e(sample)
local obs=e(N)
local T=e(g_avg)
local N=e(N_g)
sum year if e(sample)
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & e(sample)
local month_start = r(min)
sum month if year==`year_end' & e(sample)
local month_end = r(max)

* cumulative sum of gamma coefficients starting at zero in period -1
nlcom ///
	(-(_b[F1.PBD]+_b[F2.PBD])) ///
	(-_b[F1.PBD]) ///
	(0) ///
	(_b[PBD]) ///
	(_b[PBD]+_b[L1.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]+_b[L3.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]+_b[L3.PBD]+_b[L4.PBD]) ///
	, post level(95)

estimates store beta_dl_fe_crisis

* plot beta coefficients
coefplot (beta_dl_fe_crisis, ciopts(recast(rcap)) recast(connected) level(95)) ///
  	, ///
  	vertical ///
  	keep (_nl_*) ///
  	coeflabel(_nl_1 = -3 _nl_2 = -2 _nl_3 = -1 _nl_4 = 0 _nl_5 = 1 _nl_6 = 2 _nl_7 = 3 _nl_8 = 4) ///
  	xline(3.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medlarge)) ///
 	yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medlarge)) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medlarge)) ytitle ("Effect on log search intensity" " ", size(medlarge)) ///
	legend(off) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medlarge))

graph export "Fig_B2_right.pdf", replace

* ------------------------------------------------------------------------------
* Figure B.3
* multiple treatments of varying intensities (Schmidheiny/Siegloch case 4)
* changing effect window from -3,+4 to -3,+18
* estimation with distributed-lags in levels
* crisis sample
* ------------------------------------------------------------------------------

* descsribe treatment status
sum PBD

local models = ""
forvalues lag = 4(1)18 {

	* list of lagged regressors
	local laglist = ""
	forvalues l = 1(1)`lag' {
		local laglist = "`laglist' L`l'.PBD"
	}
	disp `laglist'
	
	* estimate distributed-lag model in levels with 4 to 2 leads and 18 lags
	xtreg ln_GJSI ///
		F2.PBD F.PBD ///
		PBD ///
		`laglist' ///
		frac_total_ui i.yearmonth ///
		if year <= 2011 ///
		, fe vce(cluster state)
	
	estimates store gamma_dl_`lag'_fe_crisis
	local obs=e(N)
	local T=e(g_avg)
	local N=e(N_g)
	sum year if e(sample)
	local year_start = r(min)
	local year_end = r(max)
	sum month if year==`year_start' & e(sample)
	local month_start = r(min)
	sum month if year==`year_end' & e(sample)
	local month_end = r(max)

	* cumulative lags
	local betalist = "_b[PBD]"
	local betaslist = "(_b[PBD])"
	forvalues l = 1(1)`lag' {
		local betalist = "(`betalist'+_b[L`l'.PBD])"
		local betaslist = "`betaslist' `betalist'"
	}
	disp `betaslist'

	* cumulative sum of gamma coefficients starting at zero in period -1
	nlcom ///
		(-(_b[F1.PBD]+_b[F2.PBD])) ///
		(-_b[F1.PBD]) ///
		(0) ///
		`betaslist' ///
		, post level(95)
	
	estimates store beta_dl_`lag'_fe_crisis
	generate sample_dl_`lag'_fe_crisis = e(sample)

	* list all estimated models up to here
	local models = "`models' (beta_dl_`lag'_fe_crisis, recast(connected) noci)"

	* plot beta coefficients for increasing number of lags (used in slides)
	coefplot ///
	`models' ///
  	, ///
  	vertical ///
  	keep (_nl_*) ///
  	rename(_nl_* = "") at(_coef, transform(@-4)) ///
	xlabel(-3(1)18) xline(-0.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medsmall)) xsize(7) ///
 	ylabel(-0.005(0.001)0.001) yscale(range(-.0051961  0.001)) yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medsmall)) ysize(4) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medsmall)) ytitle ("Effect on log search intensity" " ", size(medsmall)) ///
	legend(off) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medsmall))
	
	graph export "Fig_B3_`lag'.pdf", replace

}

* estimate distributed-lag model in levels with 4 leads
* use short sample as with 18 lags
* estimate distributed-lag model in levels with fixed effects with 2 leads and 4 lags
xtreg ln_GJSI ///
	F2.PBD F.PBD ///
	PBD ///
	L1.PBD L2.PBD L3.PBD L4.PBD ///
	frac_total_ui i.yearmonth ///
	if e(sample) ///
	, fe vce(cluster state)

estimates store gamma_dl_fe
local obs=e(N)
local T=e(g_avg)
local N=e(N_g)
sum year if e(sample)
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & e(sample)
local month_start = r(min)
sum month if year==`year_end' & e(sample)
local month_end = r(max)

* cumulative sum of gamma coefficients starting at zero in period -1
nlcom ///
	(-(_b[F1.PBD]+_b[F2.PBD])) ///
	(-_b[F1.PBD]) ///
	(0) ///
	(_b[PBD]) ///
	(_b[PBD]+_b[L1.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]+_b[L3.PBD]) ///
	(_b[PBD]+_b[L1.PBD]+_b[L2.PBD]+_b[L3.PBD]+_b[L4.PBD]) ///
	, post level(95)

estimates store beta_dl_4short_fe_crisis

* plot beta coefficients of all models with 4 to 18 lags plus 4 lags with 18 lags sample
coefplot ///
`models' /// 4 to 18 lags
(beta_dl_4short_fe_crisis, recast(connected) noci lpattern(dash)) /// 4 lags with 18 lags sample
, ///
vertical ///
keep (_nl_*) ///
rename(_nl_* = "") at(_coef, transform(@-4)) ///
xlabel(-3(1)18) xline(-0.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medsmall)) xsize(7) ///
ylabel(-0.005(0.001)0.001) yscale(range(-.0051961  0.001)) yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medsmall)) ysize(4) ///
graphregion(color(white)) ///
xtitle (" " "Months relative to reform", size(medsmall)) ytitle ("Effect on log search intensity" " ", size(medsmall)) ///
legend(off) ///
note(" " " ", size(medsmall))

graph export "Fig_B3.pdf", replace

* ------------------------------------------------------------------------------
* Figure B.4, left
* multiple treatments of varying intensities (Schmidheiny/Siegloch case 4)
* effect window from -3 to +4
* estimation with distributed-lags in levels and in first differences
* crisis sample
* ------------------------------------------------------------------------------

* describe treatment status and scaled treatment adoption
sum PBD
sum D.PBD

* estimate distributed-lag model in levels with 2 leads and 4 lags
* see Figure B.2 above
* stored in beta_dl_fe_crisis
* recall sample
xtsum ln_GJSI if sample_dl_fe_crisis
local obs=r(N)
local T=r(Tbar)
local N=r(n)
sum year if sample_dl_fe_crisis
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & sample_dl_fe_crisis
local month_start = r(min)
sum month if year==`year_end' & sample_dl_fe_crisis
local month_end = r(max)

* estimate distributed-lag model in first differences with 2 leads and 4 lags
* no constant in first differences, i.e. no linear trend in levels
reg D.ln_GJSI ///
	F2.D.PBD F.D.PBD ///
	D.PBD ///
	L.D.PBD L2.D.PBD L3.D.PBD L4.D.PBD ///
	D.frac_total_ui ibn.yearmonth ///
	if year <= 2011 ///
	, nocons vce(cluster state)
* Note that "ibn.yearmonth" keeps all time dummies whereas "i.yearmonth" would 
* have ommitted the first year-month as base category. We need to keep all time 
* dummies because we omit the constant. Alternatively, the constant could be
* included and interpreted as linear time trend.


* cumulative sum of gamma coefficients starting at zero in period -1
nlcom ///
	(-(_b[FD.PBD]+_b[F2D.PBD])) ///
	(-_b[FD.PBD]) ///
	(0) ///
	(_b[D1.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]) ///
	, post level(95)

estimates store beta_dl_fd_crisis

* plot beta coefficients (combined FE and FD)
coefplot ///
	(beta_dl_fe_crisis, offset(-0.05) ciopts(recast(rcap)) recast(connected) level(95) label(Fixed effects)) ///
	(beta_dl_fd_crisis, offset(0.05) color(lavender) ciopts(recast(rcap) color(lavender)) msymbol(triangle) recast(connected) level(95) label(First difference)) ///
	, ///
	vertical ///
  	keep (_nl_*) ///
  	coeflabel(_nl_1 = -3 _nl_2 = -2 _nl_3 = -1 _nl_4 = 0 _nl_5 = 1 _nl_6 = 2 _nl_7 = 3 _nl_8 = 4) ///
	xline(3.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medsmall)) xsize(5.5) ///
  	ylabel(-0.01(0.002)0.002) yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medsmall)) ysize(4.6) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medsmall)) ytitle ("Effect on log search intensity" " ", size(medsmall)) ///
	legend(region(lcolor(none))) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medsmall))
	
graph export "Fig_B4_left.pdf", replace

* ------------------------------------------------------------------------------
* Figure B.4, right
* multiple treatments of varying intensities (Schmidheiny/Siegloch case 4)
* effect window from -3 to +18
* estimation with distributed-lags in levels and in first differences
* crisis sample
* ------------------------------------------------------------------------------

* describe treatment status and scaled treatment adoption
sum PBD
sum D.PBD

* estimate distributed-lag model in levels with 2 leads and 18 lags
* see Figure B.3 above
* stored in beta_dl_18_fe_crisis
* recall sample
xtsum ln_GJSI if sample_dl_18_fe_crisis
local obs=r(N)
local T=r(Tbar)
local N=r(n)
sum year if sample_dl_18_fe_crisis
local year_start = r(min)
local year_end = r(max)
sum month if year==`year_start' & sample_dl_18_fe_crisis
local month_start = r(min)
sum month if year==`year_end' & sample_dl_18_fe_crisis
local month_end = r(max)

* estimate distributed-lag model in first differences with 2 leads and 18 lags
* no constant in first differences, i.e. no linear trend in levels
reg D.ln_GJSI ///
	F2.D.PBD F.D.PBD ///
	D.PBD ///
	L1.D.PBD L2.D.PBD L3.D.PBD L4.D.PBD L5.D.PBD L6.D.PBD L7.D.PBD L8.D.PBD L9.D.PBD ///
	L10.D.PBD L11.D.PBD L12.D.PBD L13.D.PBD L14.D.PBD L15.D.PBD L16.D.PBD L17.D.PBD L18.D.PBD ///
	D.frac_total_ui ibn.yearmonth ///
	if year <= 2011 ///
	, nocons vce(cluster state)
* Note that "ibn.yearmonth" keeps all time dummies whereas "i.yearmonth" would 
* have ommitted the first year-month as base category. We need to keep all time 
* dummies because we omit the constant. Alternatively, the constant could be
* included and interpreted as linear time trend.

estimates store gamma_dl_18_fd_crisis

* cumulative sum of gamma coefficients starting at zero in period -1
nlcom ///
	(-(_b[FD.PBD]+_b[F2D.PBD])) ///
	(-_b[FD.PBD]) ///
	(0) ///
	(_b[D1.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]+_b[L11D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]+_b[L11D.PBD]+_b[L12D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]+_b[L11D.PBD]+_b[L12D.PBD]+_b[L13D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]+_b[L11D.PBD]+_b[L12D.PBD]+_b[L13D.PBD]+_b[L14D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]+_b[L11D.PBD]+_b[L12D.PBD]+_b[L13D.PBD]+_b[L14D.PBD]+_b[L15D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]+_b[L11D.PBD]+_b[L12D.PBD]+_b[L13D.PBD]+_b[L14D.PBD]+_b[L15D.PBD]+_b[L16D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]+_b[L11D.PBD]+_b[L12D.PBD]+_b[L13D.PBD]+_b[L14D.PBD]+_b[L15D.PBD]+_b[L16D.PBD]+_b[L17D.PBD]) ///
	(_b[D1.PBD]+_b[LD.PBD]+_b[L2D.PBD]+_b[L3D.PBD]+_b[L4D.PBD]+_b[L5D.PBD]+_b[L6D.PBD]+_b[L7D.PBD]+_b[L8D.PBD]+_b[L9D.PBD]+_b[L10D.PBD]+_b[L11D.PBD]+_b[L12D.PBD]+_b[L13D.PBD]+_b[L14D.PBD]+_b[L15D.PBD]+_b[L16D.PBD]+_b[L17D.PBD]+_b[L18D.PBD]) ///
	, post level(95)

estimates store beta_dl_18_fd_crisis

* plot beta coefficients (combined FE and FD)
coefplot ///
	(beta_dl_18_fe_crisis, offset(-0.08) ciopts(recast(rcap)) recast(connected) level(95) label(Fixed effects)) ///
	(beta_dl_18_fd_crisis, offset(0.08) color(lavender) ciopts(recast(rcap) color(lavender)) msymbol(triangle) recast(connected) level(95) label(First difference)) ///
	, ///
	vertical ///
  	keep (_nl_*) rename(_nl_* = "") at(_coef, transform(@-4)) ///
	xlabel(-3(1)18) xline(-0.5, lpattern(dash) lwidth(thin) lcolor(black)) xlabel(,labsize(medsmall)) xsize(5.5) ///
  	ylabel(-0.01(0.002)0.002) yline(0, lcolor(black) lwidth(thin)) ylabel(,labsize(medsmall)) ysize(4.6) ///
 	graphregion(color(white)) ///
  	xtitle (" " "Months relative to reform", size(medsmall)) ytitle ("Effect on log search intensity" " ", size(medsmall)) ///
	legend(region(lcolor(none))) ///
	note(" " "Observations: `obs', states: `N', periods: `T' (`month_start'/`year_start' - `month_end'/`year_end').", size(medsmall))
	
graph export "Fig_B4_right.pdf", replace
